function F = cal_integral(const, q_indicator, ExportProbQ, ImportProbQ0, ImportProbQ1)

    density_omega = const.density_omega;
    E_zmQ = const.E_zmQ;
    E_zvQ = const.E_zvQ;
    E_zxQ = const.E_zxQ;

%% Calculate integrals

    integ.E_zm00 = density_omega'*(E_zmQ.*(1-ExportProbQ).*(1-ImportProbQ0).*q_indicator);
    integ.E_zm01 = density_omega'*(E_zmQ.*(1-ExportProbQ).*ImportProbQ0.*q_indicator);
    integ.E_zm10 = density_omega'*(E_zmQ.*ExportProbQ.*(1-ImportProbQ1).*q_indicator);
    integ.E_zm11 = density_omega'*(E_zmQ.*ExportProbQ.*ImportProbQ1.*q_indicator);

    integ.E_zv00 = density_omega'*(E_zvQ.*(1-ExportProbQ).*(1-ImportProbQ0).*q_indicator);
    integ.E_zv01 = density_omega'*(E_zvQ.*(1-ExportProbQ).*ImportProbQ0.*q_indicator);
    integ.E_zv10 = density_omega'*(E_zvQ.*ExportProbQ.*(1-ImportProbQ1).*q_indicator);
    integ.E_zv11 = density_omega'*(E_zvQ.*ExportProbQ.*ImportProbQ1.*q_indicator);

    integ.E_zx00 = density_omega'*(E_zxQ.*(1-ExportProbQ).*(1-ImportProbQ0).*q_indicator);
    integ.E_zx01 = density_omega'*(E_zxQ.*(1-ExportProbQ).*ImportProbQ0.*q_indicator);
    integ.E_zx10 = density_omega'*(E_zxQ.*ExportProbQ.*(1-ImportProbQ1).*q_indicator);    
    integ.E_zx11 = density_omega'*(E_zxQ.*ExportProbQ.*ImportProbQ1.*q_indicator);    
   
    
%% Return

    F = integ;

end